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Abstract 

We study spin tunneling in magnetic molecules, with special reference to 
Fes- The article aims to give a pedagogical discussion of what is meant by the 
tunneling of a spin, and how tunneling amplitudes or energy level splittings 
may be calculated using path integral and discrete phase integral methods. In 
the case of Fes , an issue of great interest is the oscillatory tunnel splittings as 
a function of applied magnetic field that have recently been observed. These 
oscillations are due to the occurrence of diabolical points in the magnetic 
field space. It is shown how this effect comes about in both the path-integral 
and the discrete phase integral methods. In the former it arises due to the 
presence of a Berry-like phase in the action, which gives rise to an interference 
between tunneling trajectories. In the latter, it arises due to the presence of 
further neighbor terms in the recursion relation for the energy eigenfunction. 
These terms give rise to turning points which have no analog in the one- 
dimensional continuum quasiclassical method. Explicit results are obtained 
for the location of the diabolical points in Fes. 

I. INTRODUCTION 

Tunneling is a basic way in which the difference between quantum and classical mechanics 
manifests itself, and even though the simplest examples of tunneling were studied right after 
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the birth of quantum mechanics, there are many other aspects of tunneling that are still the 
subject of active research today. It is generally quite difficult to calculate tunnel splittings in 
many particle systems, and even for one particle, the step from one to two spatial dimensions 
presents significant challenges and surprises UJ^j. 

In this article, I shall study the tunneling of a single spin degree of freedom. This is yet 
another instance where the statement of the problem is very simple, yet the earliest study 
of which I am aware dates to 1978 |§, a half-century after the founding of modern quantum 
theory. 

From the point of view of this volume, the greatest interest in spin tunneling lies in 
the possibility of observing meso or macroscopic quantum phenomena (MQP) in magnetic 
particles and related systems, as first proposed in the late 1980's HH. The range of activity 
over the next half-decade is well represented in a workshop proceedings || . After preliminary 
investigation, small magnetic particles appear to be attractive candidates for MQP, since 
for diameters ~ 50 A, and typical anisotropy constants, the tunneling rates appear to be 
moderately large. However, in contrast to the situation that prevails in the SQUID systems 
(for which, see the article by Han in this volume), it turns out that the classical dynamics of 
the net magnetic moment of a small particle is not fully understood theoretically, especially 
with regard to dissipation. On the experimental side, the characterization and reproducible 
fabrication of very small particles is still rather difficult, and at present even the measurement 
and modeling of the classical dynamics is not a fully solved problem. Wernsdorfer's article 
in this volume gives a good sense of the issues involved. It is this author's opinion (which 
is not necessarily shared widely), that unless the classical dynamics is fully understood, 
theoretical analyses are not likely to be pertinent, and the interpretation of experiments will 
be uncertain. 

In the last few years, however, a much more fruitful avenue for the spin tunneling has 
opened up in the area of large magnetic molecules. Several dozen molecules are presently 
under study, but the two that have yielded the most fruitful results are Mn 12 and Fe 8 . In this 
author's view, which is again not necessarily held by the majority of workers in the field, the 
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hysteresis phenomena seen in Mni2 to date fTJHl] ( c -f- review by Friedman in this volume) 
are not due to spin tunneling alone, but reflect a much richer and more complex many- 



body effect originating in the spin-phonon interaction flO||LT|j . In the case of Fes, however, 
there is unambiguous evidence of spin tunneling in the form of an oscillatory magnetic field 
dependence of the tunnel splitting, as seen by Wernsdorfer and Sessoli |T2]]. This dependence 
is highly systematic, and not easy to obtain by accident, so as an experimental signature it is 
very robust. Theoretically, these oscillations represent a very interesting difference between 
spin and massive particles. This difference can be attributed almost tautologically to the 
difference in the commutation relations. A much more visual representation of the difference 
is provided by the difference in path integrals between these systems. The spin path integral 
contains a kinetic term which has the properties of a Berry phase, which can give rise 
to interference between different spin trajectories |13|JT^| . The oscillations are a result of 
such an interference effect, and were in fact predicted theoretically without knowing of the 



relevance of the work to Fe§ ||15|| . Very briefly, here is how the effect arises. The classical 
state of a spin is defined by giving its orientation n, and the paths lie on the surface of a unit 
sphere. The sum over paths is dominated by least action paths, or instantons. For certain 
field orientations, one finds that there are two such least action paths that wind around H 
in opposite directions (Fig. 1). The real part Sr of the Euclidean action is equal for both 
paths, but the imaginary parts differ, giving rise to a relative phase equal to the Berry phase 
for the closed loop formed by the two paths. This Berry-phase is proportional to the area 
Q of the loop. Thus, the splitting A is given by 

A oc exp(-S R ) cos$, (1.1) 

where $ = JQ/2, with J being the magnitude of the spin. As the field is increased, the 
minima between which the instantons run move toward each other, and the area Q shrinks. 
Whenever <3> passes through an odd multiple of tt/2, A vanishes. In Fig. 1, we also show the 
result of a direct numerical computation of A as a function of H for the model Hamiltonian 
used in Ref. |15|, showing that the effect is real. 
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FIG. 1. Interfering instanton trajectories, and numerically computed —10 
ting for a model for Fes with H||x. 
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However, a proper discussion of spin path integrals requires the development of a fair 
amount of calculational machinery, which is not part of the standard repertoire of most 
physicists despite having been around for almost 25 years [p!6j|T7| . Indeed, in discussions 
with non-experts this author has found that several more elementary issues need to be 
clarified first. Further, path integrals are just one way to calculate tunnel splittings. In 
analogy with massive particles, there is also a phase integral or WKB method |T8|-f23f that 
can be applied to spin. For some calculational purposes, this is in fact superior to the path 
integral approach. Unfortunately, this method is also not standard textbook material, and 
while it involves far more elementary mathematics than the path integrals do, again, a fair 
amount of machinery must be developed before it can be used efficiently. 

In this article, I shall (among other things) try and give a simple treatment of both 
the path integral and discrete phase integral methods. In order to orient this discussion, 
however, it is first necessary to go back and ask what we mean by spin tunneling to begin 
with. To this end, let us consider the toy Hamiltonian 

n toy = hJ 2 x + k 2 Jl (i.2) 

with k± > k 2 > 0. J is a dimensionless spin operator with the usual commutation rules: 

[Jj, Jj] = ieijk-Jk- (1-3) 
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The first step is to understand the classical version of this problem. What is the configuration 
space, the phase space, how are the classical dynamics to be defined, and what are the 
classical states between which tunneling will take place when quantum mechanics is turned 
on? To answer these questions, we first note that a classical angular momentum (which we 
also denote by J) obeys the Poisson brackets 

{ Jj, Jj}pB = tijkJk- (1-4) 

Thus we can regard Eq. (|1.2|) as a classical Hamiltonian (with the dimensions of J and the 
fc's suitably readjusted), which defines the dynamics of the vector J through the Poisson 
brackets (p~4|) : 

dJ r -r „ ,-, . &H ,„ . 

_ = {J, W}pB = -Jx-. (1.5) 

In fact, it is clear that ( |1.5| ) is a general prescription applicable to any Hamiltonian that is 
a function of the components of J only. An immediate consequence of Eq. (|1.5|) is that 



| J . J =-2J.(jxf)=0, (1.6) 

so I J I is a constant of the motion, and configurations are completely specified by giving the 
orientation J = J/|J|, and the configuration space is the unit sphere. At the same time, 
once J is specified at any given time t — 0, Eq. (|1.5|) allows us to find it uniquely at any 
later time t > 0, as it is a first order differential equation for J. Thus the unit sphere is 
a carrier manifold in the language of modern classical mechanics, and it is also the phase 
space. In other words, for spin, configuration and phase space are one and the same. 

If we interpret Eq. ( |1.2| ) as a classical Hamiltonian, then the classical energy minima 
arise when J = ±Jz. These minima are degenerate. When the problem is made quantum 
mechanical, these two classical states will appear as two energy eigenstates with a small 
tunneling induced splitting. The problem can be approached in two stages. First of all, 
classical states with definite values of J go over into quantum mechanical states in which 
J automatically has a spread since the components of J do not commute. The states with 



J = ±Jz will acquire a zero point spread (Jj) ~ (Jy) ~ J. (One can easily find these states 
by solving and quantizing the equations of motion in the harmonic approximation. This is 
completely equivalent to the Holstein-Primakoff procedure.) In the second stage, tunneling 
mixes these states. 

With this preamble, we can pose the basic question. What is the tunnel splitting for 
a Hamiltonian such as (|1.2| ) which has two (or more) degenerate classical minima? In the 
course of answering this question, other questions arise fairly quickly. For example, from 
a purely theoretical or mathematical perspective, what are the associated wave functions? 
What are the splittings between excited pairs of levels? The answers to these can be sought 
at varying levels of rigour and quantitative accuracy. From the perspective of making contact 
with experiments, one may want to know the influence of various perturbations and "dirt 
effects". Here, an important distinction needs to be made between static and dynamic 
perturbations. For static perturbations, the problem reduces to finding the matrix elements 
of various operators between the tunneling states. The number of important operators 
in any system is generally not large, so this type of question can be moved over into the 
theory column in some sense. For time dependent perturbations, two further subclasses 
need to be recognized. If the time dependence is of the c-number type, i.e. due to a time- 
varying external field, the problem can be reduced to either a standard NMR type, or to 
a Landau-Zener-Stiickelberg type. If the perturbation is due to other dynamical degrees of 
freedom, however, the problem is much harder, and is in fact conceptually the same as that 
in investigations of dissipation and decoherence in MQP. The range of possibilities here is 
quite large in general, but in molecular magnetic systems one has the advantage of knowing 
the relevant environmental degrees of freedom with a high level of confidence, which greatly 
aids in theoretical modeling. 

In this article, we shall largely be concerned with the theoretical aspects of the problem. 
For specificity, we will focus on the Fes system, but the methodology is completely general. 
In Sec. II, we will review the salient features of the Fe 8 system, and discuss the results of 
the Wernsdorfer and Sessoli experiment, with emphasis on the oscillatory field dependence 
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and vanishing of the tunnel splitting. We will compare these results with numerical studies 
of the simplest model Hamiltonian for Fe 8 , consider the points where the splittings vanish in 
light of general quantum mechanical theorems about degeneracy, and see that these points 
form a very rich pattern in the magnetic field space. 

In Sec. Ill, we will turn to the spin-coherent-state path integral approach. These path 
integrals are much more delicate than the Feynman integral for a massive particle, and the 
mathematical subtleties of the semi-classical limit are still being researched. We will sidestep 
these points, and concentrate on the Berry phase and the quenching condition for the special 



case where the magnetic field is along the hard magnetic axis of the molecule ||15|| . This is 
the simplest case, and corresponds to a symmetric double well problem. When the field is 
not along the hard or easy exes, the problem does not have any symmetry, and instanton 
calculations, though possible, are quite involved. More quantitative calculations can be done 
with comparatively greater ease using a discrete phase integral (or WKB) method. We will 
discuss the basic idea behind this method, and then show how this method can be used to 
find tunnel splittings for Fe 8 for all orientations of magnetic field pf-E? . 



II. THE Fe 8 SYSTEM 
A. Summary of experimental facts and spin model 

The molecule Fe 8 (proper chemical formula: [Fe80 2 (OH) 12 (tacn) 6 ] 8+ ) is magnetic, and 




FIG. 2. The Fe 8 molecule. 
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forms good single crystals. It has an approximate D2 symmetry (see Fig. 2). In its low- 
est state, it is found to have a total spin of 10, arising from competing antiferromagnetic 
interactions between the Fe 3+ ions within a molecule. Spin-orbit and spin-spin interactions 
destroy complete rotational invariance, and give rise to an anisotropy with respect to the 
crystal lattice directions. A variety of experimental techniques (electron spin resonance, ac 
susceptibility, magnetic relaxation, Mossbauer spectroscopy, neutron scattering) indicates 
|28| [32| that a single molecule can be described by the Hamiltonian 

Ho = hJ 2 x + k 2 J 2 y - gf i B 3 • H, (2.1) 

with J = 10, k\ ~ 0.33 K, and hi ~ 0.22 K. In addition, there are much weaker higher order 
anisotropies. The leading fourth order anisotropy correction is 

^ = -k A {Jl + Ji), (2.2) 



with fc 4 « 2.9 x 10~ 5 K ||12j| . The anisotropy energy is equivalent to a field of ~ 2.5 T. The 
g factor is very close to 2. 

In writing the Hamiltonian for the eight coupled Fe 3+ spins in terms of a single total spin 
J as we have done, the chief assumption is that other spin multiplets are well separated in 
energy from the ground multiplet. It is very hard to do first principles calculations of the 
intramolecular, interionic exchange parameters, and even when one can do this, it is very 
hard to diagonalize the resulting Heisenberg exchange Hamiltonian. EPR experiments do 
not show any evidence of other multiplets, and so while a definite number is not known, it 
is not unreasonable to guess that other multiplets will be separated from the J = 10 ground 
multiplet by at least tens of Kelvin. Especially for the tunneling, the other multiplets can 
be safely ignored. 

Let us also briefly discuss environmental degrees of freedom which have been left out 
of the description Q2.ip , and their interaction with the spin. First, different molecules may 
interact with each other. However, the Feg molecule is very large, and in the solid, so is the 
intermolecular separation. The primary interaction between molecules is dipolar (there is 
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no intermolecular exchange, in particular), and the dipolar field on any molecule is about 
100 Oe [Q, much weaker than the intramolecular anisotropy field. Second, there is a 
spin-phonon interaction, which as in the case of Mni2, is responsible for the moderate to 



high temperature magnetization relaxation |3(| , and the dramatic hysteresis loop steps |3T 



seen in this molecule too. At low temperatures, however, very few phonons are present in 
equilibrium, and to a first approximation they may be neglected. This is especially so for 
the tunneling phenomena we will consider in this article. Third, the atomic spins couple to 
the nuclear spins-the hyperfine interaction. The biggest such coupling comes from magnetic 
nuclei in the magnetic species. For iron, the only isotope with a magnetic nucleus is 57 Fe, 
with a natural abundance of 2-2.5%. Thus, about 75% of the molecules have no nuclear spin 
in the magnetic species at all, and for those that do, the hyperfine field seen by the electronic 
spin is about 10-100 Oe. If we average the hyperfine field, and lump this along with the 
dipolar field in the form of an inhomogeneous field, we get a distribution with a width of 



about 200 Oe , far smaller than the anisotropy field. (Alternatively, we could say that the 
energy scales associated with dipolar plus hyperfine interactions and the magnetic anisotropy 
are 0.2 K and 25 K, respectively.) This approximation omits the dynamical aspects of the 
nuclear spins, which are expected to give rise to a small unquenching of the spin tunneling, 
i.e., make the transition probability nonzero. Thus, this aspect of the problem is potentially 
important for a detailed understanding of the Wernsdorfer-Sessoli data, as is the distribution 
of dipolar and hyperfine fields. 

In the rest of this article, we shall only study the pure spin problem described by Eq. ( p.l|) 
and (|2.2|) . Further, for the most part, we shall ignore the fourth order correction ( |2.2|) . This 
term is important in that even though it is a small correction to the energy of any state, 
it significantly modifies the location of the points in the magnetic field space where the 
splitting vanishes. For the conceptual problem of understanding why we get a vanishing 
splitting in the first place, however, it is sufficient to study only Ho- We shall discuss the 
effects of including Ti A is Sec. II D briefly. 
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B. Tunneling states and the Landau-Zener-Stuckelberg process 



If H = 0, Ho is exactly the toy model of Sec. I, and there are two degenerate energy 
minima at ±z. If H 7^ 0, but still in the xy plane, the classical minima move off the ±z 
directions but continue to be degenerate. One way to understand the tunneling between 
these directions is to rewrite Eq. (27T) after subtracting out a constant term k 2 J ■ J = 
k 2 J(J + 1), and explicitly putting H z = 0. This yields 



H 



-k 2 J 2 z + (h - k 2 )J 2 x - gfi B (JxH x + JyH y ). 



(2.3) 



Let us now regard the last three terms as perturbations that give rise to transitions between 
various Zeeman levels or eigenstates of J z . As usual, we denote the J z value by m. The J 2 
term gives rise to Am = 2 transitions, and thus mixes m = —10 with m = +10 via the —8, 
—6, . . ., +8 states. The H x and H y terms give rise to Am = 1 transitions, and mix m = —10 
with +10 via all intermediate levels —9 to +9 (see Fig. 3a). This picture allows us to think 
of spin tunneling in direct analogy with a particle tunneling through an energy barrier. It is 
further obvious that if we also apply a field along the z direction so as to tune the energies 
of the —10 and +9 states to resonance, we can also think of tunneling between these states. 
More generally, we can consider tunneling between a state with m = m; on the negative m 
side and m = rrif on the positive m side. 

2 r 
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m = 10 (or 9, 8, etc.) 
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(a) (b) 
FIG. 3. (a) Zeeman levels of Fes, showing Am = 1 (solid lines) and Am = 2 (dashed lines) 



transitions, (b) The Landau-Zener-Stiickelberg process. 
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At this point, it is appropriate to ask if the tunnel splitting is big enough to measure 
experimentally. Although one of the main goals of this article is to understand such split- 
tings analytically, let us temporarily cheat and diagonalize the 21 x 21 Hamiltonian matrix 
numerically with the experimental numbers for k\ and k%. When this is done, we find that 
the splitting, A ~ 10~ 8 K. We now recall that the bias or departure from degeneracy be- 
tween the m = ±10 states caused by the dipolar and hyperfine fields, which we denote by 
e, is about 0.1 K, which is enormous compared to A. Hence, left to itself, a spin will have 
almost no chance of tunneling at all. To see this, consider a two level Hamiltonian 



( e A\ 
V -A -e, 



(2.4) 



It is now easy to verify that if we start at t = in the +e level, the probability of finding 
the system in the — e level at a later time oscillates with a frequency, 2(A 2 + e 2 ) 1 / 2 , and an 
amplitude A 2 /(A 2 + e 2 ). If e 3> A, the transition probability is always very small. The 
same general conclusions apply to any two levels m; and rrif, not just the ±10 states. 

Wernsdorfer and Sessoli solve this difficulty by making use of the Landau-Zener- 



Stiickelberg (LZS) mechanism to induce tunneling transitions pi ]. A dc longitudinal field 
H z is applied so as to bring and rrif into approximate resonance. They then apply an 
additional ac longitudinal magnetic field H z (t) in the form of a triangular wave. Denoting 
the amplitude of this wave by Hq and the time period by r, we have H = \dH z /dt\ = 4Hq/t. 
Now, as H z changes with time, the energies of the —10 and +10 levels will move in opposite 
directions, and at some point in the cycle they will cross. (See Fig. 3b.) This gives rise to 
what is known as the LZS process. Basically, in the vicinity of the crossing, the energy bias 
goes to zero, and there is an appreciable chance for the spin to tunnel. The probability for 
a transition during one crossing is given by ^5 



P = l-e~\ (2.5) 

where 
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and Am = \rrif — rrii\ is the change in m in the transition. The quantity 7 is known as the 



7 < 1, and P m 7. We can further assume that because the stray fields are fluctuating, 
the phase of the molecule will be randomized and uncorrelated between successive crossings. 
If H is large enough to overcome the range of dipolar biases (but not so large as to make 
more than one Zeeman level on the positive side cross the level on the negative side), every 
molecule in the sample will undergo a crossing at some point in the cycle. Since there are 
2/t crossings per unit time, we obtain a transition probability per unit time for mj <-> m/ 
transitions, given by 



Wernsdorfer and Sessoli first saturate the sample in a large longitudinal field. This field is 
then removed and the ac field is applied, inducing LZS transitions, and causing a relaxation 
of the magnetization of the sample. By measuring the rate of this relaxation, one can obtain 
Tlzs, from which one can in turn infer A using Eq. Q2.7D and the experimental value of Hq. 
An important check on the consistency of the LZS interpretation implied by Eq. (|2.7|) is that 
the relaxation rate should be independent of the sweeping rate H . Experimentally, this is 
found to be true for H ranging from 1 mT/s to 1 T/s. 

C. Oscillatory tunnel splittings, the von Neumann- Wigner theorem, and diabolical 



It is apparent that the LZS measurements can be carried out in the presence of a trans- 
verse dc field (in the xy plane) H±, so that A can be measured as a function of H±. Naively, 
we expect that since increasing H± decreases both the energy barrier and the angle through 
which the spin must tunnel, A will increase monotonically with H±. What is actually seen 



adiabaticity parameter, and in the limit where the crossing is passed rapidly, i.e., H is large, 




(2.7) 



points 
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FIG. 4. Measured splittings [12] for Fe8 for (a) —10 «-> 10 transitions for various orientations 
of H in the xy plane, and (b) for H||x between the states m = —10 and m = 10 — n. 

experimentally is rather different (see Fig. 4). When H||y (0 = 90°), the behavior is indeed 
monotonic, but when H||x, one finds that A oscillates with H x . 

We have already stated in Sec. I that this oscillation can be understood in terms of a 
Berry phase in the spin path integral. Since the vanishing of A means that two energy levels 
of the system are exactly degenerate, it is useful to examine this result in general quantum 
mechanical terms, and from the perspective of rigorous results about when such degeneracies 
can and cannot occur. Before we turn to this, however, it is useful to look at the results of a 
numerical diagonalization of the Hamiltonian ( |2.1|) . These numerical data reveal a number 
of other properties, some of which are special to the form fl2.1|) , but others are general. In 
Fig. 5, we show the results of numerical calculation of the energies as a function of H x , for 




FIG. 5. Spectrum of the Hamiltonian (2.1) for J = 3, as a function of H x /H c , with 
H c = 2kiJ/gfJ,B- H z /H c = 0, 0.07454, and 0.1491 in (a), (b), and (c), respectively. The small ovals 
indicate narrowly avoided anticrossings that appear to be crossings on the resolution of this figure. 
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J = 3, for three different values of H z . In all three cases, H y = 0. In part (a), H z = 0, 
and problem is like a symmetric double well. We see that the lowest two energy level curves 
cross a number of times. These crossings correspond to the curve marked n — 0, i.e. to 
— 10 <-> +10 transitions in Fig. 4b. In addition, Fig. 5a also shows a number of crossings of 
higher energy levels. These crossings are difficult to see directly but their presence has been 
inferred indirectly by studying the temperature dependence of the shape of the minima in 
H x dependence of the transition rate [j36 |. 

In Fig. 5b, H z has a specific non-zero value, chosen so that (ignoring tunneling), the 
first excited state in the deeper well is degenerate with the lowest state in the shallower 
well. The problem is no longer symmetric, and one of the classical minima is lower than 
the other. Correspondingly, we see that the lowest quantum mechanical state is always non 
degenerate. However, we see from the figure that the second and third energy levels cross a 
number of times. These crossings correspond to the curve marked n = 1, i.e. to —10 <-> +9 
transitions in Fig. 4b. As seen in the experiments, the crossings in Fig. 5b are shifted by 
half a period from those in Fig. 5a. Further, as in part (a), we see crossings between yet 
higher energy levels (the fourth and fifth, e.g.). 

This pattern continues as H z is increased still further (Fig. 5c). Now the lowest two 
levels in the deeper well are nondegenerate, and the lowest crossings are between levels 3 
and 4. Compared to Fig. 5b, these crossings are shifted by yet another half-period, just as 
seen experimentally. Again, there are crossings between higher pairs of levels. 

Let us now recall the conditions under which energy levels of a quantum mechanical 
system may intersect under variation of a parameter. This is governed by the von Neumann- 



Wigner theorem [37j, which states that as a single parameter in a Hamiltonian is varied, an 
intersection of two levels is infinitely unlikely, and that level repulsion is the rule. It is useful 
to review the argument behind this theorem. Let the energies of levels in question be Ei 
and E2, which we suppose to be far from all other levels. Under an incremental perturbation 
V, the secular matrix is 



14 



'tfi + Vii v 12 ^ 

) 

V21 £?2 + ^22 ; 

with V21 = Vi2- The difference between the eigenvalues of this matrix is given by 



(2.8) 



[(£1 -E 2 + V n - K 22 ) 2 + 4|^ 12 | 2 ] 1/2 , (2.9) 

which vanishes only if 

E x + = E 2 + V 22 , V 12 = V* 2 = 0. (2. 10) 

Hence, for a general Hermitean matrix, three conditions must be satisfied for a degener- 
acy, which in general requires at least three tunable parameters. If the matrix is real and 



symmetic, the number of conditions and tunable parameters is reduced to two [38]. Degen 



eracies of the latter type, i.e. those obtained by tuning more than one parameter are known 



as diabolical, or in older terminology, as conical intersections \ 39,40]. The reason for this 
terminology is that if we denote the experimentally controllable parameters by x and y, and 
define these to be zero at the intersection, then, in its vicinity, we may expand the various 
contributions to Eq. ( |2.10| ) as 

Ei + V u - E 2 + V 22 « a x x + a y y, 
V12 ~ b x x + b y y, 

and the energy surface is given by 

E = constant ± [{a x x + a y y) 2 + (b x x + byy) 2 } 1 ^ 2 . (2-H) 

This is an elliptic double cone in the xy plane, resembling in shape an Italian toy called the 
diavolo. 

An exception to the no-crossing rule occurs when the Hamiltonian has some symmetry, 
when levels transforming differently under this symmetry can intersect. In the Fe§ moel, 
the Hamiltonian is invariant under 180° rotations about H when H||x or H||z, so there 
is such a symmetry, and states which are even or odd under the relevant operation can 
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cross 



|4T| . When H x and H z are both non-zero, however, there is no geometrical symmetry, 
and the crossings in Fig. 5(b) and (c), and the corresponding minima seen by Wernsdorfer 
and Sessoli, are non trivial instances of diabolical points. Note that the conditions of the 
theorem cited above are met no matter how one counts the parameters. If we regard the 
system as having two parameters, H x and H z , then the Hamiltonian can be chosen to be real 
by using the standard representation of the angular momentum matrices. All three angular 
momentum matrices J x , J y , and J z , cannot be made real simultaneously, however, so if we 
regard the parameter space as three dimensional (H x , H y , H z ), the Hamiltonian is complex. 
In either case, a degeneracy can only occur at an isolated point. Thus, viewed either in the 
larger H x -H z plane, or in the full three-dimensional space of magnetic fields H, all points 
of degeneracy, including those on the H x and H z cLXGS, £1X6 diabolical. 

It should be noted that in real physical systems, even when more than parameter can be 
varied, diabolical points are quite rare. Fes is remarkable in having such a rich pattern of 
intersections. Although, as already mentioned, the locations of the diabolical points depend 
sensitively on the presence of the fourth order term H4, it is of considerable theoretical 
interest that for the simpler model where the Hamiltonian is taken to be just 7{q, a number 
of results can be proved exactly P2| . (Surprisingly, the semiclassical analyses flr5| , pi| -|2"7|] 
seem to capture these results exactly at leading order in 1/J.) The first of these is for the 
location of the diabolical points. The point where the £'th level in the negative J z well (with 
£' = being the lowest level) and the £"th level in the positive one are degenerate, is at 



H y = 0, and 



H z (£',£") V\(£"-£'\ (212) 



H r 2J 

'j-n-\{£! + £" + !)] , (2.13) 



H x (£',£") 



H c J 

with n = 0, 1, . . . , 2 J - (£' + £" + 1). Here, A = k 2 /h, and H c = 2k 1 J/gfi B - Thus, the 
diabolical points lie on a perfect centered rectangular lattice in the H x -H z plane (Fig. 6). 
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FIG. 6. Diabolical points for the Hamiltonian 2.1 for J = 7/2. Many of the points are mul- 
tiply diabolical, i.e. correspond to more than one pair of simultaneously degenerate levels. This 
multiplicity is the same for all points on a rhombus and is as shown. 



Secondly, many of the points are multiply diabolical, i.e., more than one pair of levels is 
simultaneously degenerate. It is easily shown that the multiplicity is as indicated in Fig. 6: 
If we arrange the points into concentric rhombi, those on the outermost rhombus are singly 
diabolical (i.e., there is only one pair of degenerate states), those on the next rhombus are 
doubly diabolical (two pairs of degenerate states), and so on. 

These facts hint very strongly at the presence of a higher dynamical symmetry, i.e., an 
additional conserved quantity. This symmetry has not been found so far, but knowing it 
would be a tremendous advance even for real Fes, as it would be only weakly broken. For the 
same reason, it would be very useful to know the exact wavefunctions at the diabolical points, 
since that would enable one to study corrections and perturbations more systematically. 



D. Influence of higher order anisotropy perturbations 

It is of some interest to ask what happens to the diabolical points when the fourth-order 
perturbation (|2.2|) is included. Let us first investigate this question without making use of the 
specific form of the perturbation, from the general point of view of enlarging the parameter 
space of the Hamiltonian. Keeping H y = 0, we think of our Hamiltonian as depending 
on three parameters, H x , H z , and k^. The general argument about the codimension of a 
degeneracy [3|| implies that in the three-dimensional (H x , H z , k$) space, a diabolical point 
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FIG. 7. Trajectories of diabolical points under the influence of a perturbation £4. The values 
±7r and are the Berry phases associated with the adjacent contours. 

in the initial (H X ,H Z ) plane turns into a line (see Fig. 7). In the figure, we show the two 
kinds of possible behavior that are permitted by this theorem, and that are topologically 
allowed. The first, as in the line marked 'a', shows a diabolical point that continues on 
indefinitely. The second possible behavior is shown in the line marked £ b'. What appear to 
be distinct diabolical points in the k^ = plane, lie, in fact, on the same diabolical line in 
the three-dimensional space. More generally, diabolical lines can formed closed loops, but 
cannot terminate abruptly. 



A second way of viewing this matter is provided by Berry's phase [I3|. Suppose 
that at some value of fc 4 , two states \ijj a {k^ H)) and \ipb(k^ H)) are degenerate at H = 
(H x o, 0, H z q) = Ho- Let C be a small closed contour in the H x -H z plane encircling the point 
H . Berry's phase is given by 

7 (C) =z£(^ a (A;4,H)|VH^a(A:4,H))-rfH, (2.14) 

where H is now the two-dimensional vector (H X ,H Z ). As shown by Berry, 7(C) = ±7r if 
C encloses a true diabolical point, and 7(C) = if the two states merely approach each 
other very closely without ever being degenerate. [Actually, since our Hamiltonian is real 
and the parameter space (H X ,H Z ) is two-dimensional, we really only need the weaker and 
older result due to Herzberg and Longuet-Higgins PH , which states that the states \ip a ) and 
\ipb) change sign upon encircling the degeneracy: e 17 ^ = — 1. This sign change test is an 



efficient way of searching for diabolical points numerically.] Now suppose that any of the 
three parameters is changed slightly. Since the perturbation corresponding to this change 
is non-singular, |"0a(&4, H)) is a smooth function of k^ and H. It follows that if fc 4 varies 
continuously, the integrand of Eq. ( |2.14 ) can not change discontinuously. Hence, for small 
enough <5/c 4 , the phase 7(C) must continue to be what it was for fc 4 = 0, +7r, say, implying 
that C continues to encircle a degeneracy if it did so at A; 4 = 0. 

From this point of view, the behavior 'b' in Fig. 7 can only arise if 7(C) has opposite 
values for the contours encircling the two diabolical points at low values of k\. (Naturally, 
both contours must have the same sense.) The Berry phase for a contour C2 encircling both 
points is then 0, and it is then possible that for fc 4 exceeding some value k%, we can shrink 
C2 to zero, without encountering any singularity. It is obvious that this can happen only if 
the two diabolical points annihilate each other at k\. Pictorially, we can imagine "slipping" 
the contour Ci off the diabolical line by moving it above the hairpin bend in the figure. 

In Fig. 8, we show the results of a numerical calculation (performed by E. Kegecioglu), 
for J = 10, using the experimental values of k\ and k<i pertinent to Fe§. The figure shows 
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FIG. 8. Degeneracy fields for lowest two energy levels with H||x as a function of fourth order 
anisotropy. The quantity C is identical to /c 4 , and is given in Kelvin. 
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how the diabolical points on the line H z = that correspond to the degeneracy of the lowest 
two energy levels, move under the influence of the fourth order term (C = k&). We see three 
examples of a hairpin bend where two diabolical points annihilate. There are two points 
worthy of comment. First, the points that were on the surface H z = when k± = 0, continue 
to be on the same surface when fc 4 ^ 0. This can be seen as a consequence of symmetry. 
When H y = H z = 0, Tio + Tii continues to be invariant under a 180° rotation about x. Thus, 
under a change in k&, levels which cross at k& = because they have different signs under 
this symmetry operation, can continue to cross only if we continue to have H y = H z = 0. 
Second, by the time we get up to = 3 x 10 -5 K, only four diabolical points survive on 
the H x axis, and the spacing between them is nearly 50% greater than the period at = 0. 
Both these facts are in agreement with the experimental data. In fact, it was by using this 
argument in reverse, i.e., by fitting to the observed spacing that Wernsdorfer and Sessoli 
deduced the value of quoted in Sec. II A. Since the pattern and location of diabolical 
points is sensitively influenced by k 4 , this would appear to be a more reliable method of 
finding higher order anisotropy coefficients than direct EPR spectroscopy. 



III. INSTANTON CALCULATION OF TUNNEL SPLITTINGS 

In this section, we will turn to the instanton method f44jj45j of calculating the ground 
state tunnel splitting, which is based on spin coherent state path integrals. Our main goal is 
to understand the quenching effect fairly rapidly, so we will skip lightly over the subtleties 
in the spin path integral and the semiclassical limit, which are far more vexing than those 
for massive particles. (We will, however, briefly describe what these subtleties pertain to in 
subsection E.) 

A. What to calculate: the imaginary time propagator 

We use the overcomplete basis of spin coherent states {|n)} to describe our spin. The 
state |n) has maximal spin projection along the unit vector n: 
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J-n|n) = J|n). (3.1) 

We now consider the imaginary time transition amplitude 

U 21 = (nzle-^lfii), (3.2) 

where n 12 are the minima of the classical energy 

E cl (n) = (h\H\h). (3.3) 

Let us now denote the two lowest energy eigenstates by \ip±), and their energies by E±. 
where 

E± = E w ± §A, (3.4) 

with i? av being the average, and A the splitting as defined earlier. We expect that since 
these states are tunnel split states, they can be well approximated as linear combinations of 
states \ipi t 2) that are well localized around the corresponding directions n 1>2 • In other words, 

fe} = i(|^ife)), (3.5) 

with 

= 0^0, (3.6) 

(n#,)~0, (i^j). (3.7) 

Note that we can always choose the phases of so that Eq. (|3.5|) is correct as written. 

Next, we expand the amplitude U 2 i in terms of the complete set of energy eigenstates. 
As T — > oo, only the lowest two states will contribute, and we get 

E+T i /- | „/, \/„/. i~ \-E-T 



U 21 « (n 2 |^ + )(^ + |n 1 )e- ii + i + (n 2 |^_)(^|fi 1 )e 



1 * —E+T 1 * — E T 

-^a 2 a x e + — -^a 2 a x e 



a 2 ale-^ vl sinh(AT). (3.8) 
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B. The spin-coherent-state path integral 



Our goal is to calculate U 2 i via path integrals, and compare with Eq. ( j3.8|) in order 
to obtain A. The essential factor that one seeks to capture is sinh(AT). As in the case 
of massive particles we write the amplitude as a sum over paths on the unit sphere, 
weighting each path by the exponential of the action for that path: 

U 21 = r[dh]e- s ^\ (3.9) 



The paths n(r) all run from Hi at — T/2 to n 2 at T/2, and S^rifY)] is the imaginary time 
or Euclidean action for the path. (This is why the exponent is — S rather than iS/h.) The 
novel aspects for spin lie in the nature of this action, which is given by 



/■T/2 

S[n(r)} = UA[n(r)} + / E cl [n(r)]dr. (3.10) 

J-T/2 



Instead of deriving this result, we shall show that it is correct by checking that its variation 
leads to the classical equation of motion. The term «A[n(r)] is the kinetic term, and has the 
mathematical structure of a Berry phase. (The same term is responsible for the Haldane gap 
in one-dimensional antiferromagnets |47|j48[| .) One of the key properties of such a term is that 
it can not be made manifestly gauge invariant, and this fact has led to misstatements in the 
literature in which a representation arising from a particular gauge choice, and the coordinate 
singularities of spherical polar coordinates are said to be responsible for its topological 
properties. In order to avoid these pitfalls, we write it as 



A[h(r)] = i dQ, (3.11) 

J C=ti(r')— nc 



>C=h(T)-h R 

by which we mean that -4.[n(r)] is the solid angle enclosed by the closed curve formed by the 
path n(r) and a reference path n# (taken backwards) running from fii to n2 (Fig. 9). The 
reference path n# is arbitrary, but must be the same for all n(r) in the path integral. Its 
choice is equivalent to fixing the gauge. However, since we have defined A as a geometrical 
quantity, an area, it does not depend on how we choose coordinates on the unit sphere, and 
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FIG. 9. The kinetic term in the spin action, and its variation 



it is obviously nonsingular. 

The actual evaluation of A for any given closed path is often more simply done by 
using Stokes's theorem to transform it to a line integral. Using notation borrowed from 
electromagnetism, let us write dfl = B • hds, where B(n) = n, and ds is an area element. 
The line integral is / A • dh, with B = V x A. Since B = n, it is the magnetic field of a 
monopole. It is known that if we try and represent a monopole field in terms of a vector 
potential A, then A must have a singularity somewhere. If this singularity is concentrated 
into a Dirac string at the south pole, we can write 



where the curve C is regarded as being parametrized by 0. This formula is correct as long 
as C does not pass through the south pole, and provided one increments or decrements by 
2n every time one crosses the date line. 

Let us now obtain the classical equations of motion by varying the action. Consider the 
kinetic term first. Suppose we vary the path n(r). The variation 5 A is given by the area 
of the thin sliver on the sphere enclosed between the curves n(r) to n(r) + 8h(r) (Fig. 9). 
The part of this area due to the segments between r and r + 5r is given by 




(3.12) 



A(SA) 



Sh(r) x [n(r + At) - n(r)] • n(r) 




(3.13) 
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Adding up the contributions from all the segments, we find the total change 

5A = [ 5h(r) • ( ^ x n(r) ] dr. (3.14) 



dr 

The variation of the second term in Eq. ( |3.10| ) is 

dE cl 
dh 



dr. (3.15) 



Thus the variation of S can be written as an integral of the form J 5h(r) • X where X is 
something depending on n and E c \. The condition for 5S to vanish is thus X = 0, or 

dh , , <9-E r i 
iJ— x n r) + — - = 0. 3.16 
ar an 

Taking the cross product of this equation with n, and making use of the fact that n • 
(dh/dr) = 0, we get 

Zj ^ = -( nX ^-J- (3 - 17) 

This equation is exactly what we would get from Eq. ( |1.5| ) with J = Jn and the Wick 
rotation t — > —it. In other words, it is the imaginary time equation for Larmor precession 
in the effective magnetic field dE c \/dn. (It is also called the Landau-Lishitz equation in 
magnetism.) 

One consequence of Eq. (|3.17p is that E c \ is conserved along the classical path. For, 



dE r i dE r ] dh 



dr dh dr 

- 1 9Ea . (n x dE A 
J dh \ dh J 

= 0. (3.18) 



[Another way to see this is that if Eq. ( p,17|) is written in terms of spherical polar coordinates, 



it will be seen to have a Hamiltonian structure with and J cos 8 as canonically conjugate 
variables, and E c \ as the Hamiltonian.] 



24 



C. How to calculate the propagator: instantons 



The path integral ( |3.2j ) is not easy to evaluate. In the T — > oo limit, however, we can use 
the approximation of steepest descents. The dominant paths, known as instantons, are just 
those that minimize the action, i.e. they are solutions to the classical equations of motion. 
The simplest such paths consist of a single transit fro hi to n 2 . If the scale over which E c \ 
varies is V, then it follows from Eq. ( |3.17| ) that the time scale for this transit is tq ~ J/V. 
For the toy Hamiltonian ( |1.2| ), e.g., this time scale is J~ 1 {k\k'i)^ 1 l 2 . Since T — * oo, it follows 
that the spin spends most of its time near the end points n 12 , and the actual transit takes 
place in a very short time interval. (Hence the name instanton.) Further, since the equation 
of motion is autonomous, i.e., does not depend on r explicitly, in the T — > oo limit, it is 
clear that a translation of the center of the instanton yields an equally good classical path. 
Once this is realized, it is not difficult to see that one can have multi-instanton solutions, 
in which the path goes between hi and n 2 several times, with the centers of the instantons 
being widely separated on the time scale r . When all the contributions to U21 are evaluated 
in this way, one finds that the n-instanton terms give a contribution proportional to T n , and 



the full series is that of a sinh |45| . The A which is obtained in this way can be written as 

A = Dexp(-S iBBt ), (3.19) 

where Si nst is the action for a single instanton path, and D is a prefactor arising from doing 
the path-integral over small fluctuations about the instanton trajectory. If more than one 
instanton exists, we must add together the corresponding contributions from all of them. 

It is important to note that since E c \ is conserved along the instantons, and since hi 
and n 2 are minima of E c \, there can not be any real path n(r) conecting them. The only 
solution is to allow n to become complex. Correspondingly, the area A must be defined 
on the complexified unit sphere. We can take E c \ to be zero for an instanton by adjusting 
the zero of energy, so the action is just iJA. The problem is thus reduced to finding the 
instantons and the corresponding area. This is extremely simple, however. Since we can 
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write ^4[n(r)] as the line integral Q3.12| ), we do not need to find the actual time dependence 
of the instanton path, and it suffices to find the orbit on the unit sphere, i.e., 9 in terms of 
0, which can be done from energy conservation alone. 



D. Application to Fe§ 



Let us illustrate this procedure for Feg, with the Hamiltonian 7i for the case where H||x. 
If we choose the polar axis to be x (not z) , and measure the azimuthal angle in the yz plane 
from y, then we can write 



E c i( 



k\J (cosO — cos#o) + sin #sin 0. 



(3.20) 



We have defined cos^q = H/H c , with H c = 2k±J / ' g\iBi and added a constant to E c \ so that 
it vanishes at the minima (6,4>) = (#o>0) and (6 ,tt). Thus, along the instanton, E c \ = 0. 
Writing cos^o = uq, the solution of this equation gives 

M + iX 1/2 sin 0(1 -Uq-X sin 2 0) 1/2 



cos 9 



1 - A sin' 



(3.21) 



It is clear from symmetry that there are two instanton paths, which wind about x in 
opposite directions (see Fig. 1 again). We take 0(— oo) = and 0(oo) = ±7r for these paths. 
If we denote the two paths by A and B, then the real parts of their actions (= iJA) are 
equal: 



1 — A sin <p 



J 



In 



u 



In 



(1-m 2 )(1-A)+u \/A n 



(3.22) 



^l-u 2 -V\) v/T^A \^(1- U 2)(1-A)-WA, 
The imaginary parts, on the other hand are necessarily unequal, since by the interpretation 
of A as an area, 



Sb — Sa = iJ x Q, 



(3.23) 



where Q [See Eq. fl3.ll )] is the area enclosed between A and B. From Eq. (|3.12 ) we obtain, 
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Since A oc cos(JT2/2), we conclude that it vanishes whenever 

(3.25) 

in exact accord with Eq. (|2.13|) (put £ = £' = 0). Although the precise location of the 
quenching points will change if the Hamiltonian is varied, the effect is clearly general. 
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E. Tunneling prefactors 

The above discussion does not explain how the prefactor D is to be calculated. In fact, 
this calculation is somewhat more subtle for spin than it is for massive particles. If one 
evaluates the Gaussian fluctuations that yield the prefactor naively, directly following the 
massive particle case, the result is then not asymptotically correct as J — > oo |49|| . This 
point is perhaps not of great concern for numerical estimates of tunneling rates in a genuine 
physical setting, but it is nevertheless an annoying gap in the formalism. Although there do 
exist other path integral aproaches which find the splitting correctly ]5Ti| , |5"T| ] , the calculations 



are very intricate, and the simplicity seen in the massive particle case is lost. 

To describe these subtleties, let us first consider the analog of the propagator (|3.2j ) for 
a massive non relativistic particle with a geometrical position coordinate q. This is just the 
amplitude to go from a state \qi) at t — ti to a state \qj) at t = tf. Further, it is useful 
to let qi and qj be completely general, i.e., not necessarily minima of the potential energy, 
and also to consider the propagator for real time. If we write this amplitude as a Feynman 



path integral |46j , in the semi-classical limit it is again dominated by the classical paths, for 
which the action S c \ is least. Paths far away from the classical one have phases that vary 
extremely rapidly under small changes of the path, and thus these paths end up canceling 
each other. A great deal of quantum mechanics can be understood just by stringing a phase 
factor exp(iS c \/h) on the classical trajectories. If one wants to get the actual magnitudes of 
transition amplitudes, however, one must go a little further, and evaluate the integral over 
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the small fluctuations about the classical path. This gives rise to the so-called fluctuation 
or van Vleck determinant, which is a prefactor multiplying the phase factor exp(iS c \/h), c.f. 
Ref. ||52| . (The phase and prefactor correspond approximately to the eikonal and transport 
equations in the standard WKB method.) 

For spin path integrals, the fluctuation determinant is more difficult, because in contrast 
to the particle states for which (q\q') = if q 7^ q', spin coherent states are not orthogonal: 
(n|n') 7^ in general even if n ^ n'. At first sight, it seems that one should include 
discontinuous paths among the fluctuations. It turns out however, that this is not so, and 
by discrete time dissection of the path integral, the analog of the van Vleck determinant can 
be found provided one pays careful attention to the boundary conditions on the spin path. 



This has been done by several authors ||53|-|56|j , but the results do not seem to be widely 
known. It seems to this author, that once the van Vleck prefactor is understood, tunneling 
prefactors should be calculable with the same degree of ease as for massive particles |45| , 



but except for some work in Ref. |55| this does not seem to have been widely appreciated 
yet. 



IV. THE DISCRETE PHASE INTEGRAL METHOD 

A. Basic formalism: local Bloch waves 

The basic idea of the DPI method is to try and solve Schrodinger's equation in the J z 
basis as a recursion relation or difference equation. For constant coefficients, linear difference 
and differential equations can both be solved in terms of exponentials. For slowly varying 
coefficients, the WKB approach is a powerful one for differential equations, and many of 
the ideas used there can be carried over to difference equations. With this in mind, let us 
suppose |^) is an eigenstate of H with energy E. Then, with J z \m) = m\m), {m\ip) = C m , 
{m\TL\m) = w m , and (m\H\m') = t m ^ m i (m ^ m'), we have 

^ ^ ^m,nC?i ^mPm EC m , (4.1) 
n 
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where the prime on the sum indicates that the term n = m is to be omitted. 

We can think of Eq. ( |4.1| ) as a tight binding model for an electron in a one-dimensional 
lattice with sites labelled by m, and slowly varying on-site energies (w m ), nearest-neighbor 
(t m m±1 ) hopping terms, next-nearest-neighbor (t mjm ±2) hopping terms, and so on. In most 
interesting problems, hopping to very distant neighbors is negligible, and the recursion 
relation involves only a handful of terms. Since we can think of dynamics in this model 
in terms of wavepackets, it is clear that there is a generalization of the usual continuum 
quasiclassical or phase integral method to the lattice case. This is the DPI method. 

More specifically, the DPI method is applicable to a recursion relation such as ( |4.1|) 
whenever the w m and t mim ± a vary sufficiently slowly with m. If these quantities were in- 
dependent of m, the solutions to Eq. (|4. 1| ) would be Bloch waves C m = exp(zgm), with an 
energy 

E = w m + 2t rn ^ l+ i cosg + 2t m>m+2 cos2g = E(q), (4.2) 

To save writing, it is sometimes useful to identify w m = t m ^ mi and we shall use the notations 
w m and t mim interchangably. If for fixed a, the t mim+a vary slowly with m (where the 
meaning of this term remains to be made precise), we expect it to be a good approximation 
to introduce a local Bloch wavevector, q(m), and write C m as an exponential e**, whose 
phase $ accumulates approximately as the integral of q(m) with increasing m, in exactly 
the same way that in the continuum quasiclassical method in one dimension, one writes the 
wavefunction as exp(iS(x)/h), and approximates S(x) as the integral of the local momentum 
p(x). 



Previous work with the DPI method 18.19,|57|,p0|,pl| has been limited to the case where 



the recursion relation has only three terms, i.e., only nearest neighbor hopping is present. 
As discussed by Braun, the DPI approximation has been employed in many problems in 
quantum mechanics where the Schrodinger equation turns into a three-term recursion re- 
lation in a suitable basis. All the types of problems as in the continuum case can then be 
treated — Bohr-Sommerfeld quantization, barrier penetration, tunneling in symmetric double 
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wells, etc. In addition, one can also use the method to give asymptotic solutions for various 
recursion relations of mathematical physics, such as those for the Mathieu equation, Hermite 
polynomials, Bessel functions, and so on. The general procedures are well known and simple 
to state. For any E, one solves the Hamilton- Jacob i and transport equations to obtain q(m) 
and v(m), and writes C m as a linear combinations of the independent solutions that result. 
The interesting features all arise from a single fact — that the DPI approximation breaks 
down at the so-called turning points. These are points where v (m) vanishes. One must 
relate the DPI solutions on opposite sides of the turning point by connection formulas, and 
the solution of all the various types of problems mentioned above depends on judicious use 
of these formulas. 

In the Feg problem, the recursion relation involves five terms. The diagonal terms (w m ) 
arise from the Jf and J Z H Z parts of H, the i m>m ±i terms from the J X H X part, and the t m>m ±2 
terms from the part. This gives rise to new features over and above the three term case. 
In particular, we encounter nonclassical turning points, i.e., turning points at m values other 
than those at the limits of the classically allowed motion. It is these turning points that 
give rise to oscillatory tunnel splittings, so that this effect is absent in systems described by 
three-term recursion relations. 

The fundamental requirement for a quasiclassical approach to work is that w m and t m , m ± Q 
(a = 1, 2) vary slowly enough with m that we can find smooth continuum approximants 
w(m) and t a (m), such that whenever m is an eigenvalue of J z , we have 

w{m) = w m , (4.3) 

t a (m) = (t m , m+Q + tm,m-a)/2, « = 1, 2. (4.4) 

We further demand that 

with m/J being treated as quantity of order 1. Problems for which these conditions cannot 
be met are not amenable to the DPI method. It is not difficult to see that for Eqs. Q2.1| ) 



and (|2.2j ), these conditions will hold in the semiclassical limit J ^> 1. 
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Given these conditions, the basic approximation, which readers will recognize from the 
continuum case, is to write the wavefunction as a linear combination of the quasiclassical 
forms 

C m ~ exp (i [ q(m')drr/) , (4.6) 



'v{rn) 

where q(m) and v (m) obey the equations 

E = w(m) + 2t\(m) cosg + 2t 2 (?7i) cos(2g) = H sc (q, m), (4.7) 
v(m) = d7i sc /dq = — 2 smq(m)(ti(m) + 4t 2 (m) cosg(m)). (4-8) 

Equations ([4.7|) and (^4.8|) are the lattice analogs of the eikonal and transport equations. 
Equation (^4.6|) represents the first two terms in an expansion of log C m in powers of 1/ J. 

B. Turning points and connection formulas 

The basic DPI approximation fails whenever v(m) = &H sc (q,m)/dq = 0, because then 
it diverges. We shall call all such points turning points in analogy with the continuum case. 
In contrast to that case, however, we will find that turning points are not just the limits of 
the classical motion for a given energy, once the notion of the classically accessible region is 
suitably understood. 

Since we must also obey the eikonal equation ( [4.7D in addition to the condition v{m) = 0, 
at a turning point both m and q are determined if E is given. Setting v — in Eq. ( (4.8|) . 
we see that we must have either q = 0, or q = ir, or q = g*(m), where 

cosg*(m) = — ti(m)/4t 2 (^)- (4.9) 

Substituting these values of q in the eikonal equation, we see that a turning point arises 
whenever 

E = U (m), U„{m), or U m (m), (4.10) 

where, 
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U (m) = H sc (0,m) = w(m) + 2t 1 (m) + 2t i (m), (4.11) 

£4(m) = H sc (ir,m) = w(m) - 2*i(m) + 2t 2 (m), (4.12) 

t 2 (m) 

U*(m) = H sc {q*,m) = w{m) - 2t 2 {m) - -j^r. (4.13) 

Note that q*(m) may be complex for some m, but since cosg* is always real, U* is real 
for all m. We shall refer to these three energy curves as critical curves. We shall see that 
they collectively play the same role as the potential energy in the continuum quasiclassical 
method. 

To better understand the turning points, let us assume that t\ < 0, and t 2 > 0. [This is 
the case for the Hamiltonian (|2.1[) . We can always arrange for t\ to be negative by means 
of the gauge transformation C m — > (— l) m C m . Thus there is only one other case to be 
considered, namely, t\ < 0, t 2 < 0. This is discussed in Ref. |58|].] It then follows that 
U n > Uq, and that 

U (m) - U*{m) = (m) + 4t 2 (m)) 2 > 0. (4.14) 

4:t 2 {m) 

Secondly, let us think of T-[ sc (q,m) for fixed m as an energy band curve. Then U n is always 
the upper band edge, while the lower band edge is either Uq or [/# according as whether 
—t\/At 2 is greater than or lesser than 1. To deal with this possibility, it pays to introduce a 
dual labeling scheme for all three curves U , U n , and U*. We write U^(m) = U + (m), and 

U (m) = Ui(m), U*(m) = C/_(m), if g* e (0,tt), (4.15) 
U (m) = U-(m), U*(m) = U f {m), if g* ^ (0, tt). (4.16) 

The subscripts + and — denote upper and lower band edges, while the subscripts i and / 
denote internal and forbidden respectively, since in the first case above, U lies inside the 
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u = u_ 



-m 



u.=u 



-m o / - m 1 m 

FIG. 10. Critical energy curves for the Fes Hamiltonian when H||x, showing the dual labeling 
scheme. In the right hand figure, the region near the left minimum of U$ is magnified, showing the 
two turning points —mt and —mi. 

energy band, while in the second case, U* lies outside. As examples of these curves for a 
symmetric recursion relation, we show those for Fes in Fig. 10, along with a magnified view 
of the lower left hand portion. 

Turning points where E = U + , or E = £/_ when LL = Uq, are analogous to those 
encountered in the continuum quasiclassical method, since the energy lies at a limit of the 
classically allowed range for the value of m in question. Points where E = U- when [/_ = U* 
are physically analogous, but mathematically different since the value of q c is neither nor 
7i. Points where E = U (see the energy E\ in Fig. 10, e.g.) are novel in that the energy 
is inside the classically allowed range for m c , but the mathematical form of the connection 
formulas is identical to the case E = C7_ = Uq since q c = 0. Most interesting are the 
turning points with E = Uf (the point m = —mi in Fig. 10, for instance), since now the 
energy is outside the allowed range for m = m c , and the value of q c is therefore necessarily 
complex. These points lie "under the barrier" and turn out to be the ones of importance for 
understanding oscillatory tunnel splittings. The derivation of connection formulas at these 
turning points is not particularly difficult, but quite lengthy. In fact, it is quite lengthy 
to even quote the connection formulas themselves, and we refer readers to Ref. |58 for the 
details. The key point is that a solution with a purely imaginary q on side of the turning 
point, representing, let us say a solution that decays with growing m, turns on the other 
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side into a solution with a complex q with a nonzero real part, representing a solution that 
decays with an oscillating envelope as m increases. These kinds of wavefunctions can simply 
never arise in conventional one-dimensional continuum problems. 

From this point on, the exercise of calculating tunnel splittings is exactly as in the 
continuum case. One demands that the wavefunctions decay as m — > ±00 and applies 
connection formulas at the turning points to continue the wavefunction into the interior. 
Finally near m = 0, one demands that the wavefunctions so found from the two sides should 
agree. This gives the eigenvalue condition. The one new complication (besides the new 
turning points) is that at each m, there are four solutions to the Hamilton- Jacobi equation 
instead of two in the continuum case. Asm-> —00, all four q's are pure imaginary, but of 
these only the two which are negative imaginary can be kept. Thus at every point m, one 
must keep track of two DPI solutions. An exception to this arises in the case of symmetric 



double wells (H||x), where Herring's formula |59| simplifies matters, and one can work 
purely with a wavefunction that decays away from one of the wells in both directions. This 
formula is sufficiently important to be worth discussing separately. We do this in the next 
subsection, present results that follow from its application in the succeeding subsection, and 
for the asymmetric case in the one following that. 

C. The symmetric case; Herring's formula 

It is useful to recall a basic formula for tunnel splittings in symmetric potentials 
[V(— x) = V(x)] for massive particles |60 |. Let the minima be at x = ±a, and let ipo{x) be 



a wavefunction localized entirely in the right hand well, with an energy E . This wavefunc- 
tion obviously does not obey Schrodinger's equation near the left hand well, and we could 
imagine modifying the potential suitably in that region so that the wavefunction continues 
to decay wtih decreasing x even in the vicinity of the left well. Since we will never need 
ipo{~ a )i the precise way in which this is done is not important. 

If Eq is well below the barrier, then we expect there to be two states ip s and ipa, with 
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energies E\ and E 2 , both very close to E , and with wavefunctions that are very accurately 
given by 

The product ip (x)ip (— x) is everywhere exponentially small, and so therefore, if the wave- 
function i^o(x) is normalized, so are ip_ %a . For the same reason, 

^ ^ Q {x)^ a {x)dx = -j=. (4.18) 

The differential equations obeyed by ipo and ip a in the region x > are 

- (k 2 /2m)^(x) + V(x)Mx) = EoMx), (4-19) 
-(h 2 /2m)£(x) + V{x)ij a {x) = E 2 ^ a (x). (4.20) 

We multiply the first equation by ip a , the second by ip , subtract, and integrate from to 
oo. This yields 

h 2 /-oo. 



h roo 

E 2 -E = V2— / faift - ^a)da 
2m Jo 



x 

Jo 



, 7i 

ft 2 



m 



^o(0)Vo(0), (4.21) 



where we have used the facts that ^0,0(00) = 0, V'o(O) = 0, and V'a(O) = V^dW- A similar 
calculation yields E\ — E = — {E 2 — E ), so that the energy splitting, A is given by 

2k 2 

A = —^(0)^(0). (4.22) 

This is Herring's formula. It has a nice and physically appealing interpretation in terms of 
the probability current at x — 0. 

If we now use the WKB method to find V'o(O) and ^0(0)5 keeping in mind the normal- 
ization, for the nth pair of levels, we get 



A n = g n — exp 

7T 



-o' h 



(4.23) 
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Here, u is the small oscillation frequency in the well, [V'^aj/m] 1 / 2 , ±a' n are the classical 
turning points for the nth energy level pair with the mean value E w (n + ^)hu, and 



9* 



nl 



(4.24) 



Note that if g n were unity, we would have the formula of Ref. [BO]. The corrections are small: 
go = (n/e) 1 / 2 ~ 1.075, g\ ~ 1.028, gi ~ 1.017, and so on. Stirling's formula for n\ shows 
that g n — > 1 as n — > oo. These corrections are well known to many workers, but not known 
to many others, and they arise from the fact that as n becomes smaller, the use of a linear 
turning point formula at the turning point is increasingly inaccurate. The correct procedure 
is to use quadratic turning point formulas [^]. A pedagogical discussion of this matter may 



be found in Ref. Iplj . 

One can repeat the above argument for the discrete case making the obvious modifications 



as needed. The analogue of Eq. ( 4.22 ), Herring's formula, is 



A 



^o,iCo(Ci — C-i) + £0,2^0 (C2 — O-2) + t-i t i(Cf — C^) 



2 t 



Cl-C 2 , 



4 t_s 1 [CiCs 

2 ' 2 V 2 2 



C_lC_3 1 . 

2 2/ 



integer J, 
half-integer J. 



(4.25) 



When it is combined with the DPI approximation, one obtains for both integer and half- 
integer j pufl , 

a u g n 



2tt 



exp 1 



q(m')dmj + c.c. 



(4.26) 



The similarity to the continuous case is striking. The only point worth remarking is the 
possibility of having the Bloch vector q(m) be complex in part or all of the tunneling region. 
In Eq. ( |4.26j ), we have written the formula for the Fes problem, where there are two g's with 
positive imaginary parts and equal and opposite real parts in the tunneling region. 

It may be useful at this point to make a general remark about how tunneling prefactors 
are obtained in the DPI method. Equations ( |4.23| ) and ( |4.26| ) contain these prefactors 



naturally, and it is apparent that they arise as a consequence of the l/yv(m) normalization 
in the general DPI form and the connection formulas at turning points. Provided one can 
solve the transport equation and apply the connection formulas correctly, these factors are 
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relatively easy to obtain. (We shall see that this continues to be true in the asymmetric 
case.) Historically, however, one error which has been made is to use linear turning point 
formulas all throughout, as in Ref. ||60|| . As already mentioned, this leads one to miss the 
curvature factors g n . This error is widely repeated, perhaps because of the reputation of 
this text, and perhaps because it is numerically insignificant. From a logical point of view, 
however, there is no reason not to incorporate the curvature corrections, and they must be 
kept if one wants an answer for the tunneling amplitude itself (and not just its logarithm) 
that is properly asymptotic to the true answer as h — > 0, or as J — > oo. 



D. Results for Fe§ in symmetric case [26] 

To apply the general results to Fes, we need the explicit forms for w(m), and t a (m). It 
is convenient to measure energies (including u ) in units of ki J 2 , and introduce the scaled 
variable fi = m/J, where 

J=J + \. (4.27) 

In these variables, when H||x, we get (with h x = JH X /JH C ) 

w(m) = (1 + A)(l - m 2 )/2, (4.28) 
tl {m) = -h x {l-^) l l\ (4.29) 
t 2 (m) = (l-A)(l-/i 2 )/4. (4.30) 

The actual evaluation of the integrals is mostly a matter of careful but straightforward 
analysis, and the final answer for the splitting can be written as 

A n = ^W-u;oF n +5 e - r °cosA n , (4.31) 

where, 

u, = 2J[kMl -hlo)} 1/2 , (4-32) 

^1/2(1 _ fc2\3/2 



"X 
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r = J 



In 



1 - ul + v / A N 



i-w^-Va/ Vi -a 



:ln 



l- M 2)(l-A) + ^v / A N 



(1 A) -/^VXy 



A», = max < 0, Ti J 1 



In Eqs. ( ggg[ - gg§) , A = k 2 /h, and 

JH X 

jH r ; 



nn 



h, 



H x 



Recall that H c = 2k\.J /qhb- 

It is an immediate consequence of these results that A n vanishes when 



H c J 



J 



(4.34) 
(4.35) 



(4.36) 



(4.37) 



with I = n, n + 1, . . ., 2 J — n — 1. These are precisely the results quoted in Sec. II. 



E. Results for Fe§ in asymmetric case [27] 

When H has components both along x and z, the critical curves are no longer symmetric 
(Fig. 11), and the problem cannot be reduced to just quoting a splitting. Suppose we consider 
an energy E, as drawn in Fig. 11, and suppose that it leads to turning points at m' a , m' b , m' c 
on the left hand side, and m" < m'l < m" on the right. (We denote quantities pertaining to 
the left hand solution or the left hand side of the well by either a single prime or a suffix 




m 



-m* U* 

FIG. 11. Critical energy curves for the Fes Hamiltonian when H has both x and z components. 
In the right hand figure, the region near the left minimum of Uq is magnified, showing the various 
turning points. 
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— , and corresponding right-hand-well quantities by a double prime or a + suffix.) A wave- 
function C' m which decays as m — > — oo will possess the following additional characteristics. 
It will either oscillate or continue to have the same exponential character in the classically 
allowed region m' a < m < m' b . In the region just to the right of m' b it will consist of a decay- 
ing part and a growing part. The new feature will be encountered at m! c where E = U*. For 
m > m' c , both the growing and decaying parts will acquire oscillatory envelopes. Similar 
remarks apply to the right side wavefunction C^. The eigenvalue condition will be obtained 
by matching the wavefunction in the central region m' c <C m <C m". 

Let us now denote the left well minimum by m -, the small oscillation frequency in that 
well by uj -, and the value of U (m -) by E_. We write 

E = E_ + + i)w _. (4.38) 

In the same way we define analogous quantities for the right hand well. We expect that 
tunneling will be significant only when the energy E is such that it matches a level in both 
wells, say, level numbers £' and £", where £ = is the ground level. Since the problem 
is harmonic near the two potential wells, we can write the energy levels using a harmonic 
oscillator formula, and with this in mind, we define an offset e', 

u' = £' + — , (4.39) 
No- 
where £' is a positive integer, and e' lies in the interval (—1/2, 1/2)uj ~- e" is similarly defined. 

Our remarks above mean that mixing between the two wells will be significant only when 

e' and e" are both small. This is indeed the case, and the matching or eigenvalue condition 

turns out to be 

e'e" = ^[A{£',£")} 2 , (4.40) 

where A(£', £") is an exponentially small tunneling amplitude that we shall give shortly. The 
essential point emerges if we define 

e = \{e' + e") = E - \(e_ + E + + {£' + \)u> _ + {£" + ±V +) , (4.41) 
5 = e"-e'= (E_ -E + + {£' + \)uj^ - {£" + |)cu 0+ ) . (4.42) 
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Then, the eigenvalue condition Eq. (|4.40|) can be rewritten as 



e= ±\[5 2 + A 2 (f,f')] 1/2 - (4-43) 



^TLS — ^ 



(4.44) 



These are of course, the eigenvalues of the two-level-system Hamiltonian 

1 5 A(f,T)^ 

v A(f,r) -5 J ' 

exactly as we should expect. The quantity e is the energy measured from a convenient 
reference point, while S, which depends on the fields H x , H z , and the quantum numbers £' 
and £" of the states whose mixing is being examined, is the bias or offset between these energy 
levels in the absence of tunneling. Equation ( |4.45| ) below gives the tunneling amplitude 



between these levels when the bias is small, i.e., when the two levels are in approximate 
degeneracy. Note that although this amplitude is defined even for relatively large biases — 
biases comparable to the intrawell spacings u Q ± — and indeed is not very sensitive to the 
value of the bias, the concept of tunneling is physically sensible and useful only when the 
bias is comparable to or less than the amplitude A. If 5 ^> A, we get e ~ ±5/2, i.e., e" ~ 5, 
e' A 2 / 5, or the other way around. Choosing the first way, we find that E « E_ + (£'+^)ujq^, 
and a wavefunction essentially given by (1 A 2 /5) T , i.e., localized in the left well, with 
negligible mixing with the right well. 

It remains to give the expression for A(i', £"). This is, 

A(f ,£") = -(^^) 1/2 (^o-^o+) 1/2 e- rG cosA c , (4.45) 

7T 

with Tq being the total Gamow factor 

r m b 

Tq = \Im qi(m)\ dm, (4.46) 

J m' 

b 

and A c being the phase integral 

A c = / \Req(m)\dm. (4.47) 
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In Eq. (|4.46|) the subscript 1 on q(m) means that we are to take that solution of the Hamilton- 
Jacobi equation which goes to zero at m' b and m'^. In Eq. ( |4.47| ) , on the other hand, we can 
choose any of the four solutions for the wavevector, since they are of the form 

q(m) = ±m(m) ± x{ m )i (4.48) 

where /t(m) and x{ m ) are both real, and the signs can be chosen independently from the 
two ± options. 

The application of these results to Fes again requires doing a certain number of integrals. 
The problem of greatest interest is the location of the diabolical points, and for that we need 
only solve the conditions 5 = A = 0. The latter condition can only come about when A c is 
an odd multiple of 7r/2, so in fact we need only give formulas for S and A c . We find 

5(h z , £',£") = 4M, + - £") - ^(£' + £" + l) Cl (h x ) + 0(J- 3 ), (4.49) 



where 



l-h 2 x + \(l + 2hp 
Cl (h x ) = , (4.50) 



with h z = JH Z /JH C , and /io = (1 — h 2 x ) 1 / 2 . For A c , we have up to order J°, 



A c = - 
2 



2 J - (£' + £") - 2J- 



(4.51) 



If we ignore the correction Ci(h x ), the conditions 5 = A = are thus precisely those quoted 
in Sec. II. 
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